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Experiments with shadowgraph visualization reveal a subcritical transition to a hexagonal convection pattern in thin liquid 
layers that have a free upper surface and are heated from below. The measured critical Marangoni number (84) and observation 
of hysteresis (3%) agree with theory. In some experiments, imperfect bifurcation is observed and is attributed to deterministic 
forcing caused in part by the lateral boundaries in the experiment. 
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The onset of motion in heated fluid layers with a free 
upper surface has eluded complete understanding ever 
since Benard’s investigation Q of these flows established 
thermal convection as a paradigm for pattern formation 
in nonequilibrium systems [^. Rayleigh’s analysis 
of this problem assumed that buoyancy effects, which 
are always present in layers heated from below, caused 
convection, but the threshold that Rayleigh predicted 
did not agree with Benard’s observations. Forty years 
elapsed before it was recognized that the instability ob¬ 
served in Benard’s studies was caused not by buoyancy 
but by surface tension gradients as characterized by 
the Marangoni number M (see Fig. |^). Linear theory 
yields onset at Mc=80. Weakly nonlinear theory [||J^] 
predicts a subcritical (hysteretic) transition to a hexago¬ 
nal pattern. Only a single experimental investigation ||] 
has systematically examined the onset of convection for 
layers sufficiently thin such that surface tension forces 
dominate over buoyancy. That experiment revealed a 
primary transition to a concentric roll pattern at values 
of M that decreased as the fluid layers became thinner; 
for the thinnest layers studied, rolls emerged at M an 
order of magnitude smaller than Me from theory. 

In this Letter we present evidence for a well-defined 
primary transition in surface-tension-driven Benard 
(Marangoni) convection experiments designed so that 
surface tension forces dominate over buoyancy to a 
greater extent than in previous investigations. We ob¬ 
serve a hysteretic bifurcation to a defect-free array of 
hexagonal cells; this bifurcation is modeled by an ampli¬ 
tude equation, which permits comparison to both linear 
and weakly nonlinear stability theory. We also observe 
hexagons to arise from an imperfect bifurcation where the 
hysteresis disappears; this bifurcation is described quali¬ 
tatively with the addition of a deterministic forcing term 
to the amplitude equation. In our experiments, surface 
tension effects are 40 times larger than buoyancy effects, 
ie., M/i?«40, where the Rayleigh number is defined as 
R=gP A Td? IVK, with liquid expansion coefficient (3 and 


gravitational acceleration g. A necessary condition for 
the flow to be surface-tension-dominated is M/R>1 [|j; 
previous experiments attained ||^. 



FIG. 1. Cross section of our cylindrical convection cell. 
The dimensionless control parameter is the Marangoni num¬ 
ber M=aT ATd/pvK, where aT=\d(j/dT\, and o, p, v, k are, 
respectively, the liquid surface tension, density, kinematic vis¬ 
cosity and thermal diffusivity. The assumption of conductive 
heat transport is used to obtain the mean temperature across 
the liquid layer A T={Tb — Tt) /(I + H~^) with the Biot num¬ 
ber H=kgd/kdg defined in terms of the thermal conductivities 
k and kg of the liquid and gas, respectively. 

The experiments are performed on a purified silicone 
oil layer (d=0.0419±0.0005 cm) that is bounded from be¬ 
low by I-cm-thick gold-coated aluminum mirror (Fig. |^) 
[p^ . A uniform air layer (dg=0.0455±0.0008 cm) lies be¬ 
tween the oil layer and a 1-mm-thick sapphire window. 
The oil is confined by a Teflon sidewall ring of inner di¬ 
ameter 4.53±0.01 cm. Because nonuniform wetting can 
cause large relative changes in the liquid layer thickness 
near the sidewall, a polyethersulfone annular buffer of 
thickness 50 gm suppresses convection adjacent to the 
sidewall jTl| ] . The inner diameter of the buffer determines 
the radius to height ratio r=45.6±0.1 of the convecting 
region. The oil and the air layer depths each vary by less 
than 1% over the central 70% of the convecting region, 
as measured both mechanically and interferometrically. 
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Use of a purified silicone oil (96.7% hexacosamethyl- 
dodecasiloxane) avoids both condensation and cross- 
diffusive effects 1^^ that can affect pattern formation. 
A temperature gradient is imposed by water-cooling the 
window to a temperature T(=13.320±0.005°C and by 
computer-controlled-heating the mirror to a temperature 
Tb that fluctuates less than ±0.0005°C. For sufficiently 
small Tf, — Tf, the surface tension a{T) at the liquid-gas 
interface is uniform; however, with Tb — Tt sufficiently 
large, instability causes surface tension variations that 
drive flow in the bulk. The shadowgraph technique is 
used to detect onset and to visualize patterns. Images 
are digitized and background subtracted to improve the 
signal-to-noise ratio. The time scale in the experiment is 
set by the vertical diffusion time, ty=d'^/ k=1.9 s. 




FIG. 2. The abrupt onset of hexagons in Marangoni con¬ 
vection. (a) Just prior to onset, weak convection rolls develop 
at the boundary for £=-5.6x10“®. (b) A hexagonal pattern 
fills the entire convection apparatus for e=-2.5xl0“®. 


vection pattern disappears. The front between hexagonal 
and conductive states can remain stationary indefinitely 
(or at least for times long compared to the horizontal 
diffusion time) for states like Fig. ^(b) or (c) if e is held 
constant. 
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FIG. 3. Return to the conductive state for decreas¬ 
ing e. Below onset hexagonal convection persists at (a) 
£=-2.60xl0“^ (b) e=-2.82xl0“^ and (c) e=-2.96xl0“^ be¬ 
fore disappearing at (d) e=-3.20x 10“^. 


Figure || demonstrates that the conductive state un¬ 
dergoes an abrupt transition to hexagons as M is in¬ 
creased slowly (dM/dt=10“‘* in units of Just prior 

to onset, weak circular convection rolls arise near the 
boundary [upper left and lower right in Fig. ||(a)]; we 
believe these rolls are driven by static forcing due to the 
slight mismatch of thermal conductivity between the an¬ 
nular buffer and the liquid. Convection cells first appear 
within a portion of the boundary rolls after an increase 
of 2 X 10“® in e={M/Mc)-l. (My is determined from the 
experiments, as will be described.) Additional hexagons 
then nucleate from the initial cells and propagate as a 
traveling front, invading the apparatus until the entire 
flow domain is filled with the hexagonal pattern [Fig. 
11(b)] . The resulting pattern is nearly free from defects 
since the lattice is grown from a single “seed crystal” at 
the boundary. The front propagates across the appara¬ 
tus in approximately 900t„, a time short compared to the 
horizontal diffusion time 4r^tj,f«8000t„. 

On decreasing e quasistatically, hexagons persist at pa¬ 
rameter values below that for the first appearance of cells 
(Fig. ^). The transition to the conductive state occurs 
gradually; the control parameter typically must be decre¬ 
mented through a range of e before the hexagonal con¬ 


Near onset, the hexagonal pattern arises from the 
interaction of three plane wave (roll) solutions whose 
wavevectors have a magnitude equal to the critical 
wavenumber and differ in angle by 2tt/3 [ p^ . The evo¬ 
lution of the pattern can then be described by a Landau 
equation for the amplitude A 

A = €A + aA'^ - A^ + f, (1) 

with a>0 and / a constant that can account for deter¬ 
ministic forcing. In some cases, the coefficients in Eq. (1) 
can be computed from the full fluid equations jl^]. The 
existence of hexagons requires thus, the bifurcation 
from the conductive state must be subcritical. The so¬ 
lutions for hexagonal convection and for conduction are 
both linearly stable over a range of parameter: £„<£< 0 
with £^=-0^/4 (the conductive state is linearly unstable 
for e>0). 

Equation (1) is a variational model that exhibits relax- 
ational time dependence governed by a potential function 
[p^ , which we first consider for /=0. Over the parameter 
range where both conduction and convection are stable, 
each state corresponds to a minimum of the potential; 
one state represents the global minimum while the other 
state, the metastable phase, represents a local minimum. 
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The potential varies as e changes; at a parameter value 
(the Maxwell point) both states have equal values of the 
potential. As e passes through Cm, the states exchange 
the roles of global stability/metastability. For Eq. (1), 
the conductive state is globally stable for e<em=|eo, and 
metastable for em<e< 0. 



FIG. 4. Hysteresis at the onset of Marangoni convection is 
demonstrated by a plot of the Fourier mode amplitude A from 
shadowgraph images vs e. Convection appears suddenly with 
slowly increasing e (+) and persists below onset for slowly 

decreasing e (triangles). A ht to the convective branch (-) 

yields Mc=83.6 (Arc=1.65°C), which we use to compute e. 


ratus. With decreasing e, hexagonal convection can be¬ 
come metastable; however, the range of parameter values 
where hexagons are metastable is nearly an order of mag¬ 
nitude smaller than the region of metastability for con¬ 
duction. This suggests that the transition back to con¬ 
duction will be more sensitive to small spatial variations 
in e due to nonuniformities in the depths of both liquid 
and gas layers. Thus, for values of e near the metastable 
region of hexagons, the front will move in stages to spread 
the conductive state across the apparatus, as shown in 
Figs. ||(b) and|^(c). 



FIG. 5. Imperfect bifurcation in Marangoni convection, (a) 

Data (x) are compared to Eq. (1) with /=1.3/c (-). The 

parabola from Fig. 4 (-) is also shown, (b) Shadowgraph 

image of weak convective flow at e=-5.30x10“^ in the pres¬ 
ence of signihcant deterministic forcing. 


To compare the experimental observations to the 
model, we compute two-dimensional spatial power spec¬ 
tra from shadowgraph images. The spectra are az- 
imuthally averaged and normalized to the variance of the 
image intensity. The mean position of the fundamental 
spectral peak yields the wavenumber 1.90±0.02 (nondi- 
mensionalized by d); linear stability analysis predicts a 
critical wavenumber of 1.99 [^. The wavenumber is in¬ 
dependent of e for the range investigated. The amplitude 
in Fig. 1^ is the square root of the power contained in the 
spectral peak at the fundamental wavenumber. 

Figure]^ demonstrates that the experimental observa¬ 
tions illustrated in Figs. ^ and ^ are consistent with Eq. 
(1). Hexagonal convection amplitudes for increasing and 
decreasing e near the bifurcation are fit by a parabola, 
as suggested by (1) with /=0; from this fit we estimate 
Mc=83.6 with a precision of ±0.5 in M. The uncer¬ 
tainty in the accuracy is ±11 in M, primarily due to the 
uncertainty in the thermal properties for the silicone oil. 
From Fig. we also estimate ea=-3.2±0.3xl0“^ and 
em=-2.8±0.3xl0“^ lll^ . For increasing e, the conduc¬ 
tive state shown in Fig. ^(a) is deep within the metastable 
regime when the initial onset occurs. The weak convec¬ 
tion roll at the boundary provides a sufficient perturba¬ 
tion to push the system over the potential barrier, and 
the front between the two states propagates to spread the 
globally stable state (hexagons) across the entire appa¬ 


In some cases, convection appears without hysteresis. 
This situation arises, for example, in experiments where 
e is repeatedly increased and decreased, causing conduc¬ 
tion and convection to alternate. As the number of cy¬ 
cles increases, hysteresis is observed at smaller values of 
e and for a smaller range of e. The convective onset oc¬ 
curs continuously [Fig. ^a)] after a sufficient number of 
cycles; this number varies from 3 to 15 between different 
experimental runs. Our observations indicate that deter¬ 
ministic forcing, which increases slowly (on a time scale 
much longer than the horizontal diffusion time), causes 
an imperfect bifurcation [|^. Equation (1) models im¬ 
perfect bifurcation with /yfO; for f>fc=a^/27, hysteresis 
disappears at the onset of hexagons. In this regime, Eq. 
(1) qualitatively describes the amplitudes measured from 
our experiments [Fig. ||(a)]; the difference between the 
model and the data suggests that a more complex form 
for / [e.g., /(e)] is necessary for quantitative agreement. 

The physical origin of the forcing that causes imper¬ 
fect bifurcation has not been determined definitively; 
however, the observation of low amplitude rolls paral¬ 
lel to the lateral boundary [Fig. ||(b)] suggests sidewall 
boundaries are affecting the flow. Similar structures arise 
in buoyancy-driven convection with intentional thermal 
forcing at the sidewall [|^, although forcing for surface- 
tension-driven flow is probably more complex because 
some deformation of the free surface at the boundary is 
unavoidably present due to nonuniform contact line pin- 
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ning at the sidewall. The low amplitude flows become in¬ 
creasingly cellular away from the boundaries and toward 
the center of the apparatus; moreover, with increasing 
e above onset, the rolls at the sidewall are supplanted 
by hexagonal cells as the amplitudes arising from imper¬ 
fect bifurcation approach the amplitudes observed during 
hysteretic onset [Fig. ||(a)]. 

Our experimental studies of onset confirm the pre¬ 
dictions of theory and suggest an explanation for the 
puzzling disagreement between previous experiments and 
theory for the century-old problem of surface-tension- 
driven Benard convection. Our determination of Mc=84 
is in reasonable agreement with M^. from linear theory 
1^^ , and our observation of subcritical bifurcation is in 
accord with weakly nonlinear theory. Our finding of 
3.2% hysteresis sets a standard for comparison to non¬ 
linear theories, whose estimates of hysteresis range from 
0.2% 01 to 2.3% 0. Observation of imperfect bifurcation 
demonstrates the sensitivity of the primary instability in 
Marangoni convection to perturbations; the appearance 
of rolls before hexagons at M«Mc in previous experi¬ 
ments 0 may well be due to this sensitivity. 
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